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Abstract 

We study the zero-temperature phase diagram of the 2D quantum J^-^'^-J^^^ spin- 1/2 
anisotropic Heisenberg model on the square lattice. In particular, the effects of the anisotropy A on 
the z-aligned Neel and (collinear) stripe states, as well as on the xy-planar-aligned Neel and collinear 
stripe states, are examined. All four of these quasiclassical states are chosen in turn as model states 
on top of which we systematically include the quantum correlations using a coupled cluster method 
analysis carried out to very high orders. We find strong evidence for two quantum triple points 
(QTP's) at (A^ = -0.10±0.15, Jf/Ji = 0.505±0.015) and (A^ = 2.05±0.15, Jf/Ji = 0.530±0.015), 
between which an intermediate magnetically-disordered phase emerges to separate the quasiclassi- 
cal Neel and stripe collinear phases. Above the upper QTP (A > 2.0) we find a direct first-order 
phase transition between the Neel and stripe phases, exactly as for the classical case. The z-aligned 
and xy-planar-aligned phases meet precisely at A = 1, also as for the classical case. For all values 
of the anisotropy parameter between those of the two QTP's there exists a narrow range of values 
of J2/J1, a'^^(A) < J2/J1 < a'^^(A), centered near the point of maximum classical frustration, 
J2/J1 = 5, for which the intermediate phase exists. This range is widest precisely at the isotropic 
point, A = 1, where a^i(l) = 0.44±0.01 and 0:^2(1) = 0.59±0.01. The two QTP's are characterized 
by values A = A'^ at which a'^i(A^) = a^^iA"). 

PACS numbers: 75.10.Jm, 75.30.Gw, 75.40.-s, 75.50.Ee 
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I. INTRODUCTION 



The exchange interactions that lead to collective magnetic behavior are clearly of purely 
quantum-mechanical origin. Nevertheless, the underlying quantum nature has often safely 
been ignored in describing, at least at the qualitative level, many magnetic phenomena of 
interest in the past. On the other hand, the investigation of magnetic systems and magnetic 
phenomena where the intrinsically quantal effects play a dominant role, and hence have to 
be accounted for in detail, has evolved in recent years to become a burgeoning area at the 
forefront of condensed matter theory. Thus, the investigation of quantum magnets and their 
phase transitions, both quantum and thermal, has developed into an extremely active area 
of research. 

From the experimental viewpoint major impetus has come both from the discovery of 
high-temperature superconductors and, since then, from the ever-increasing ability of mate- 
rials scientists to fabricate a by now bewildering array of novel magnetic systems of reduced 
dimensionality, which display interesting quantum phenomena.- While high-temperature 
superconductivity has raised the question of the link between the mechanism of supercon- 
ductivity in the cuprates, for example, and spin fluctuations and magnetic order in one- 
dimensional (ID) and two-dimensional (2D) spin-half antiferromagnets, the new magnetic 
materials exhibit a wealth of new quantum phenomena of enormous interest in their own 
right. 

For example, in ID systems, the universal paradigm of Tomonaga-Luttinger liquicl^'^ 
behavior has occupied a key position of interest, since Fermi liquid theory breaks down 
in ID. More generally, in all restricted geometries the interplay between reduced dimen- 
sionality, competing interactions and strong quantum fluctuations, generates a plethora of 
new states of condensed matter beyond the usual states of quasiclassical long-range order 
(LRO). Thus, for high-temperature superconductivity, for example, it is suggested- that 
quantum spin fluctuation and frustration due to doping could lead to the collapse of the 
2D Neel-ordered antiferromagnetic phase present at zero doping, and that this could be the 
clue for the superconducting behavior. This, and many similar experimental observations 
for other magnetic materials of reduced dimensionality, has intensified the study of order- 
disorder quantum phase transitions. Thus, low- dimensional quantum antiferromagnets have 
attracted much recent attention as model systems in which strong quantum fluctuations 
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might be able to destroy magnetic LRO in the ground state (GS). In the present paper we 
consider a system of ^ oo spin-1/2 particles on a spatially isotropic 2D square lattice. 

The spin-1/2 Heisenberg antiferromagnet with only nearest-neighbor (NN) bonds, all of 
equal strength, exhibits magnetic LRO at zero temperature on such bipartite lattices as the 
square lattice considered here. A key mechanism that can then destroy the LRO for such 
systems, with a given lattice and spins of a given spin quantum number s, is the introduction 
of competing or frustrating bonds on top of the NN bonds. The interested reader is referred 
to Refs. [ija] for a more detailed discussion of 2D spin systems in general. 



An archetypal model of the above type that has attracted much theoretical attention in 
recent years (see, e.g., Refs. gIQ 

yyiiolQQQQQIielQQH) is the 2d spin-1/2 J1-J2 

model on a square lattice with both NN and next-nearest-neighbor (NNN) antiferromagnetic 
interactions, with strength Ji > and J2 > respectively. The NN bonds Ji > promote 
Neel antiferromagnetic order, while the NNN bonds J2 > act to frustrate or compete with 
this order. All such frustrated quantum magnets continue to be of great theoretical interest 
because of the possible spin-hquid and other such novel magnetically disordered phases 
that they can exhibit (and see, e.g., Ref. 20|]). The recent syntheses of magnetic materials 
that can be well described by the spin-1/2 J1-J2 model on the 2D square lattice, such as 
the undoped precursors to the high-temperature superconducting cuprates for small J2/J1 
values, VOM0O4 for intermediate J2/ Ji values,— and Li2VOSi04 for large J2/J1 values,— 
has fuelled further theoretical interest in the model. 

The properties of the spin-1/2 J1-J2 model on the 2D square lattice are well understood 
in the limits when J2 = or Ji = 0. For the case when J2 = 0, and the classical GS 
is perfectly Neel-ordered, the quantum fluctuations are not sufficiently strong enough to 
destroy the Neel LRO, although the staggered magnetization is reduced to about 61% of 
its classical value. Indeed, the best estimates for this order parameter are 61.4 ± 0.1% 
from quantum Monte Carlo studies,— 63.5% from exact diagonalizations of small clusters,— 
61.4 ± 0.2% from series expansions,— 61.5 ± 0.5% from the coupled cluster method (CCM) 
employed here,— i^"^ and 61.4% from third-order spin-wave theory.— Clearly, they all agree 
remarkably well in this J2 = limit. The opposite limit of large J2 is a classic example^ of the 
phenomenon of order by disorder.— 1^ Thus, in the case where Ji — with J2 7^ and fixed, 
the two sublattices each order antiferromagnetically at the classical level, but in directions 
which are independent of each other. This degeneracy is lifted by quantum fiuctuations and 
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the GS becomes magnetically ordered collinearly as a stripe phase consisting of successive 
alternating rows (or columns) of parallel spins. 

For intermediate values of J2/ Ji it is now widely accepted that the quantum spin-1/2 
J1-J2 model on the 2D square lattice has a ground-state (gs) phase diagram showing the 
above two phases with quasiclassical LRO (viz., a Neel-ordered (vr, vr) phase at smaller 
values of J2I J\i and a collinear stripe-ordered phase of the columnar (vr, 0) or row (0, vr) 
type at larger values of J^lJx)-, separated by an intermediate quantum paramagnetic phase 
without magnetic LRO in the parameter regime a^^ < J2/J1 < a'^^, where a'^^ ~ 0.4 and 
a'^^ ~ 0.6. The precise nature of the intermediate magnetically-disordered phase is still 
not fully resolved. Suggested candidates include a homogeneous spin-liquid state of various 
types with no broken symmetry (see, e.g., Ref. [3]), or a valence-bond solid (VBS) phase 
with some broken symmetry. Possible spin-liquid states include a resonating-valence-bond 
(RVB) state proposed by Anderson,- which has been supported more recently by variational 
quantum Monte Carlo studies.— Other studies^'^'^'^'^ have supported a spontaneously 
dimerized state for the intermediate phase with both translational and rotational symmetry 
broken, and thus representing a columnar VBS phase. Yet other studies^*^ have supported 
instead a plaquette VBS state for the intermediate phase, with translational symmetry 
broken but with rotational symmetry preserved. 

There has also been considerable discussion in recent years as to whether the quan- 
tum phase transition between the quasiclassical Neel phase and the magnetically disordered 
(intermediate paramagnetic) phase in the spin-1/2 J1-J2 model on the 2D square lattice 
is first-order or of continuous second-order type. A particularly intriguing suggestion by 
Senthil et al.^ is that there is a second-order phase transition in the model between the 
Neel state and the intermediate disordered state (which these authors argue is a VBS state), 
which is not described by a Ginzburg-Landau-type critical theory, but is rather described 
in terms of a deconfined quantum critical point. Such direct second-order quantum phase 
transitions between two states with different broken symmetries, and which are hence char- 
acterized by two seemingly independent order parameters, are difficult to understand within 
the standard critical theory approach of Ginzburg and Landau, as we indicate below. 

Thus, the competition between two such distinct kinds of quantum order associated with 
different broken symmetries would lead generically in the Ginzburg-Landau scenario to one of 
only three possibilities: (i) a first-order transition between the two states, (ii) an intermediate 
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region of co-existence between both phases with both kinds of order present, or (iii) a region 
of intermediate phase with neither of the orders of these two phases present. A direct 
second-order transition between states of different broken symmetries is only permissible 
within the standard Ginzburg-Landau critical theory if it arises by an accidental fine-tuning 
of the disparate order parameters to a multicritical point. Thus, for the spin-1/2 J1-J2 
model on the 2D square lattice and its quantum phase transition suggested by Senthil et 
al.,— it would require the completely accidental coincidence (or near coincidence) of the 
point where the magnetic order parameter (i.e., the staggered magnetization) vanishes for 
the Neel phase with the point where the dimer order parameter vanishes for the VBS phase. 
Since each of these phases has a different broken symmetry (viz., spin-rotation symmetry 
for the Neel phase and the lattice symmetry for the VBS phase), one would naively expect 
that each transition is described by its own independent order parameter (i.e., the staggered 
magnetization for the Neel phase and the dimer order parameter for the VBS phase) and 
that the two transitions should hence be mutually independent. 

By contrast, the "deconfined" type of quantum phase transition postulated by Senthil 
et air- permits direct second-order quantum phase transitions between such states with 
different forms of broken symmetry. In their scenario the quantum critical points still sep- 
arate phases characterized by order parameters of the conventional (i.e., in their language, 
"confining" ) kind, but their proposed new critical theory involves fractional degrees of free- 
dom (viz., spinous for the spin-1/2 J1-J2 niodel on the 2D square lattice) that interact via 
an emergent gauge field. For our specific example the order parameters of both the Neel 
and VBS phases discussed above are represented in terms of the spinous, which themselves 
become "deconfined" exactly at the critical point. The postulate that the spinous are the 
fundamental constituents of both order parameters then affords a natural explanation for 
the direct second-order phase transition between two states of the system that otherwise 
seem very different on the basis of their broken symmetries. 

We note, however, that the deconfined phase transition theory of Senthil et al^ is still 
the subject of controversy. Other authors believe that the phase transition in the spin-1/2 
J1-J2 niodel on the 2D square lattice from the Neel phase to the intermediate magnetically- 
disordered phase need not be due to a deconfinement of spinous. For example, Sirker et 
alM have argued on the basis of both spin-wave theory and numerical results from series 
expansion analyses, that this transition is more likely to be a (weakly) first-order transition 
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between the Neel phase and a VBS phase with columnar dimerization. Other authors have 
also proposed other, perhaps less radical, mechanisms to explain such second-order phase 
transitions (if they exist) and their seeming disagreement (except by accidental fine tuning) 
with Ginzburg-Landau theory. What seems clearly to be a minimal requirement is that 
the order parameters of the two phases with different broken symmetry should be related 
in some way. Thus, a Ginzburg-Landau-type theory can only be preserved if it contains 
additional terms in the effective theory that represent interactions between the two order 
parameters. For example, just such an effective theory has been proposed for the 2D spin- 
1/2 J1-J2 model on the square lattice by Sushkov et al.^ and further discussed by Sirker 
et alM 

From the classical viewpoint frustrated models often exhibit "accidental" degeneracy, and 
the degree of such degeneracy, which can vary enormously, has become widely viewed as a 
measure of the frustration. Among the effects that can act to lift any such degeneracy are 
thermal fluctuations, quantum fluctuations, and such "perturbations" as spin-orbit interac- 
tions, spin-lattice couplings, further neglected exchange terms, and impurities, all of which 
might be present in actual materials. In the present paper we focus particular attention on 
the role of quantum fluctuations. From the quantum viewpoint such frustrated quantum 
magnets as the spin-1/2 J1-J2 model on the 2D square lattice often have ground states that 
are macroscopically degenerate. This feature leads naturally to an increased sensitivity of 
the underlying Hamiltonian to the presence of small perturbations. In particular, the pres- 
ence in real systems that are well characterised by the J1-J2 model, of anisotropics, either 
in spin space or in real space, naturally raises the issue of how robust are the properties of 
the model against any such perturbations. 

Combining the above two viewpoints, it is clear that it is of particular interest in the study 
of frustrated quantum magnets to focus special attention on the mechanisms or parameters 
that are available to us to "tune" or vary the quantum fluctuations that play such a key role 
in determining their gs phase structures. Apart from changing the spin quantum number 
or the dimensionality and lattice type of the system, or tuning the relative strengths of the 
competing exchange interactions, another key mechanism is the introduction of anisotropy 
into the existing exchange bonds. Such anisotropy can be either in real space^'^'^'^'^'^ 

46 47 48 4'3 

or m spin space.— i^i^i^ 

In order to investigate the effect in real space an interesting generalization of the pure 
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J1-J2 model has been introduced recently by Nersesyan and Tsvelik— and further studied 
by other groups including ourselves.— '^'^'^'^ This generalization, the so-called Ji-J[-J2 
model, introduces a spatial anisotropy into the 2D J1-J2 model on the square lattice by 
allowing the NN bonds to have different strengths Ji and J[ in the two orthogonal spatial 
lattice dimensions, while keeping all of the NNN bonds across the diagonals to have the 
same strength J2. In previous work of our own^^^ on this Ji-J[-J2 model we studied 
the effect of the coupling J[ on the quasiclassical Neel-ordered and stripe-ordered phases 
for both the spin-1/2 and spin-1 cases. For the spin-1/2 case,— we found the surprising 
and novel result that there exists a quantum triple point below which there is a second- 
order phase transition between the quasiclassical Neel and columnar stripe-ordered phases 
with magnetic LRO, whereas only above this point are these two phases separated by the 
intermediate magnetically disordered phase seen in the pure spin-1/2 J1-J2 model on the 2D 
square lattice (i.e., at J[ = Ji). We found that the quantum critical points for both of the 
quasiclassical phases with magnetic LRO increase as the coupling ratio J'l/Ji is increased, 
and an intermediate phase with no magnetic LRO emerges only when J'l/ Ji ^ 0.6, with 
strong indications of a quantum triple point at J[/ Ji = 0.60 ±0.03, J2/J1 = 0.33 ±0.02. For 
J(/Ji = 1, the results agree with the previously known results of the J1-J2 model described 
above. 

In the present paper we generalize the spin-1/2 J1-J2 model on the 2D square lattice in 
a different direction by allowing the bonds to become anisotropic in spin space rather than 
in real space. Such spin anisotropy is relevant experimentally as well as theoretically, since 
it is likely to be present, if only weakly, in any real material. Furthermore, the intermediate 
magnetically-disordered phase is likely to be particularly sensitive to any tuning of the 
quantum fluctuations, as we have seen above in the case of spatial anisotropy. Indeed, other 
evidence indicates that the intermediate phase might even disappear altogether in certain 
situations, such as increasing the dimensionality or the spin quantum number. 

Thus, for example, the influence of frustration and quantum fluctuations on the magnetic 
ordering in the GS of the spin-1/2 J1-J2 model on the body-centered cubic (bcc) lattice has 
been studied using exact diagonalization of small lattices and linear spin-wave theory,^- and 
also by using linked-cluster series expansions. Contrary to the results for the corresponding 
model on the square lattice, it was found for the bcc lattice that frustration and quantum 
fluctuations do not lead to a quantum disordered phase for strong frustration. Rather, the 
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results of all approaches suggest a first-order quantum phase transition at a value J2/J1 ~ 
0.70 from the quasiclassical Neel phase at low J2 to a quasiclassical coUinear phase at large 
J2. Similarly, the intermediate phase can also disappear when the spin quantum number s 
is increased for the J1-J2 model on the 2D square lattice. Thus, we^^ found no evidence for 
a magnetically disordered state (for larger values of J2I J\) for the s = 1 case, by contrast 
with the s = 1/2) case.— Instead, we found a quantum tricritical point in the s = 1 case of 
the Ji-J[-J2 model on the 2D square lattice at J'JJi = 0.66 ± 0.03, J2/J1 = 0.35 ± 0.02, 
where a line of second-order phase transitions between the quasiclassical Neel and columar 
stripe-ordered phases (for J[/ Ji < 0.66) meets a line of first-order phase transitions between 
the same two phases (for J[l Ji > 0.66). 

As in our previous work— involving the effect of spatial anisotropy on the spin-1/2 
and spin-1 J1-J2 models on the 2D square lattice, we again employ the coupled cluster 
method (CCM) to investigate now the effect on the same model of spin anisotropy. The 
CCM is one of the most powerful techniques in microscopic quantum many-body theory.— 1^ 
It has been applied successfully to many quantum magnets.— '^'^'^'^'^'^ It is capable of 
calculating with high accuracy the ground- and excited-state properties of spin systems. In 
particular, it is an effective tool for studying highly frustrated quantum magnets, where such 
other numerical methods as the quantum Monte Carlo method and the exact diagonalization 
method are often severely limited in practice, e.g., by the "minus-sign problem" and the very 
small sizes of the spin systems that can be handled in practice with available computing 
resources, respectively. 

II. THE MODEL 

The usual 2D spin-1/2 J1-J2 model is an isotropic Heisenberg model on a square lattice 
with two kinds of exchange bonds, with strength Ji for the NN bonds along both the row 
and the column directions, and with strength J2 for the NNN bonds along the diagonals, as 
shown in Fig. [TJ^a). Here we generalize the model by including an anisotropy in spin space 
in both the NN and NNN bonds. We are aware of only a very few earlier investigations 
with a similar goal.— ^^"^ The two most detailed have studied the extreme limits where 
either the frustrating NNN interaction becomes anisotropic but the NN interaction remains 
isotropic^ (viz., the Ji-J^^^ model) and the opposite case where the NN interaction be- 
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FIG. 1: (a) The J^^^-J^^^ Heisenberg model; — Ji; J2; (b) and (c) z-aligned states for 

the Neel and stripe columnar phases respectively; (d) and (e) planar x-aligned states for the Neel 
and stripe columnar phases respectively. Arrows in (b), (c), (d) and (e) represent spins situated 
on the sites of the square lattice [symbolized by • in (a)]. 

comes anisotropic but the NNN interaction remains isotropic^ (viz., the J^^^-J2 model). 
In real materials one might expect both exchange interactions to become anisotropic. To 
our knowledge the only study of this case^^ (viz., the J^^^-J^^^ model) has been done 
using the rather crude tool of linear spin-wave theory (LSWT), from which it is notori- 
ously difficult to draw any firm quantitative conclusions about the positions of the gs phase 
boundaries of a system. It is equally difficult to use LSWT to predict with confidence either 
the number of phases present in the gs phase diagram or the nature of the quantum phase 
transitions between them. We comment further on the application of spin-wave theory to 
the J1-J2 model and its generalizations in Sec. |Vl The aim of the present paper is to use the 
CCM, as a much more accurate many-body tool, to investigate the spin-1/2 J^^^-J^^^ 
model on the 2D square lattice. 

In order to keep the size of the parameter space manageable the anisotropy parameter A 
is assumed to be the same in both exchange terms, thus yielding the so-called J^^^-J^^^ 
model, whose Hamiltonian is described by 



(id) 
{(i,k)) 



k) 5 



(1) 



where the sums over (i, j) and ((i, k)) run over all NN and NNN pairs respectively, counting 
each bond once and once only. We are interested only in the case of competing antiferro- 
magnetic bonds, Ji > and J2 > 0, and henceforth, for all of the results shown in Sec. 



10 



llVt we set Ji = 1. Similarly, we shall be interested essentially only in the region A > 
(although for reasons discussed below in Sec. [IV]we shall show results also for small negative 
values of A). 

This model has two types of classical antiferromagnetic ground states, namely a 2;-aligned 
state for A > 1 and an xy-planar-aligned state for < A < 1. Since all directions in the 
xy-plane in spin space are equivalent, we may choose the direction arbitrarily to be the x- 
direction, say. Both of these z- aligned and x- aligned states further divide into a Neel (vr, vr) 
state and stripe states (columnar stripe (vr, 0) and row stripe (0,7r)), the spin orientations 
of which are shown in Figs. [l]^b,c,d,e) accordingly. There is clearly a symmetry under the 
interchange of rows and columns, which implies that we need only consider the columnar 
stripe states. The (first-order) classical phase transition occurs at J2 = with the Neel 
states being the classical GS for J2 < | Ji, and the columnar stripe states being the classical 
GS for J2 > |Ji. 



III. THE COUPLED CLUSTER METHOD 



We briefly outline the CCM formalism (and see Refs. |27||52|53||y|55||5d|57|i58||59| for 
further details). The first step of any CCM calculation is to choose a normalized model 
(or reference) state 1$) which can act as a cyclic vector with respect to a complete set of 
mutually commuting multi-configurational creation operators, C/ = (C/)^- The index / 
here is a set-index that labels the many-particle configuration created in the state C/|$). 
The requirements are that any many-particle state can be written exactly and uniquely as 
a linear combination of the states {C/|$)}, together with the conditions. 



= = C7 1$) V/ ^ ; C+ = 1 



(2) 



[Ci,Cj] = 0=[CY,Cj]. (3) 
The Schrodinger equations for the many-body ground-state (gs) ket and bra states are 



(4a) 
(4b) 



respectively, with normalization chosen such that = 1 [i.e., with (\^| = ((\1/|\1/)) ^(^|], 

and with itself satisfying the intermediate normalization condition (^I^P) = 1 = (*^*|^*)- 



11 



In terms of the set {|$); C*/}, the CCM employs the exponential parametrization 

1^) =6^1$), S =J2^iCf (5a) 

for the exact gs ket energy eigenstate. Its counterpart for the exact gs bra energy eigenstate 
is chosen as 

(^| = (<l>|^e^ S =l + Zi^,SiCY. (5b) 

It is important to note that while the parametrizations of Eqs. f l5a|) and (I5bp are not man- 
ifestly Hermitian-conjugate, they do preserve the important Hellmann-Feynman theorem at 
all levels of approximation (viz., when the complete set of many-particle configurations {/} 
is truncated). Furthermore the amplitudes {Si,Si) form canonically conjugate pairs in a 
time-dependent version of the CCM, by contrast with the pairs {Si,Sj), coming from a 
manifestly Hermitian-conjugate representation for (^| = (($|e'^^e'^|<I>)^^($|e'^\ that are not 
canonically conjugate to one another.— 

The static gs CCM correlation operators, S and S, contain the real c-number correlation 
coefficients, Sj and Sj, that need to be calculated. Clearly, once the coefficients {Si,Si} 
are known, all other gs properties of the many-body system can be derived from them. To 
find the gs correlation coefficients we simply insert the parametrizations of Eqs. (ISal b) into 
the Schrodinger equations (Hal b) and project onto the complete sets of states and 
C/|$), respectively. Completely equivalently, we may simply demand that the gs energy 
expectation value, H = (\E'|if|\l/), is minimized with respect to the entire set {iS/,iS/}. In 
either case we are easily led to the equations 

(<l>|C7e-^//e^|<I)) = 0; V/ ^ , (6a) 
($|^e-^[i/,C+]e^|$) = 0; V/ ^ , (6b) 

which we then solve for the set {Sj, Sj}. Equation fl6a|) also shows that the gs energy at the 
stationary point has the simple form 

E = E{{Si}) = {^\e-^He^\^) . (7) 

It is important to realize that this (bi-)variational formulation does not necessarily lead to 
an upper bound for E when the summations for 5* and 5* in Eqs. fl5al b) are truncated. 
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due to the lack of manifest Hermiticity when such approximations are made. Nonetheless, 
one can prove^'^ that the important Hellmann-Feynman theorem is preserved in all such 
approximat ions . 

We note that Eq. fl5a|) represents a coupled set of non-linear multinomial equations for the 
c-number correlation coefficients {Si}. The nested commutator expansion of the similarity- 
transformed Hamiltonian, 

e-'He' = H+[H,S] + ^[[H, 5], 5] + ■ ■ ■ , (8) 

and the fact that all of the individual components of S in the expansion of Eq. (ISal) commute 
with one another by construction [and see Eq. ([3])], together imply that each element of S in 
Eq. (!5a|) is linked directly to the Hamiltonian in each of the terms in Eq. ([8]). Thus, each of 
the coupled equations ( l6al) is of Goldstone linked- duster type. In turn, this guarantees that 
all extensive variables, such as the energy, scale linearly with particle number, N . Thus, at 
any level of approximation obtained by truncation in the summations on the index I in Eqs. 
fl5a]b). we may (and do) always work from the outset in the limit iV — >^ oo of an infinite 
system. 

Furthermore, each of the linked-cluster equations fl6a|) is of finite length when expanded, 
since the otherwise infinite series of Eq. ([8]) will always terminate at a finite order, provided 
only (as is usually the case, including that of the Hamiltonian considered here) that each 
term in the Hamiltonian, if, contains a finite number of single-particle destruction operators 
defined with respect to the reference (vacuum) state |$). Hence the COM parametrization 
naturally leads to a workable scheme that can be computationally implemented in a very 
efficient manner. 

Before discussing the possible CCM truncation schemes, we note that it is very convenient 
to treat the spins on each lattice site in a chosen model state |$) as equivalent. In order to 
do so we introduce a different local quantization axis and a correspondingly different set of 
spin coordinates on each site, so that all spins, whatever their original orientations in |$) 
in a global spin-coordinate system, align along the negative ^-direction, say, in these local 
spin coordinates. This can always be done by defining a suitable rotation in spin space of 
the global spin coordinates at each lattice site. Such rotations are canonical transformations 
that leave the spin commutation relations unchanged. In these local spin axes where the 
configuration indices / simply become a set of lattice site indices, / ^ {fci, /c2, ■ ■ ■ fcm}, the 
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generalized multi-configurational creation operators Cj are simple products of single spin- 
raising operators, C/ ■^ki-^k2 ' ' ' ^tm' where = si ± is^. , and (s^, si, si) are the usual 
SU(2) spin operators on lattice site k. For the quasiclassical magnetically-ordered states 
that we calculate here, the order parameter is the sublattice magnetization, M, which is 
given within our local spin coordinates defined above as 

1 ^ 

k=l 

It is usually convenient to take the classical ground states as our (initial) choices for the 
model state |$). Hence, we may choose here either a Neel state or a (columnar) stripe state 
for 1$). Each of these can be further sub-divided into a z-aligned choice or a planar (say, 
x-aligned) choice, which we expect to be appropriate for A > 1 and |A| < 1 respectively 
on purely classical grounds. We present results below in Sec. HV] based on all four of these 
classical ground states as choices for 1$). 

Clearly the CCM formalism is exact when one includes all possible multi-spin configura- 
tions / in the sums in Eqs. flSaf b) for the cluster correlation operators S and S. In practice, 
however, truncations are needed. As in much of our previous work for spin-half models 
we employ here the so-called LSUBn scheme,— '^'^'^'^'^'^'^'^ in which all possible multi- 
spin-fiip correlations over different locales on the lattice defined by n or fewer contiguous 
lattice sites are retained. (Two sites are defined to be contiguous here if they are NN sites on 
the lattice.) The numbers of such fundamental configurations (viz., those that are distinct 
under the symmetries of the Hamiltonian and of the model state |$)) that are retained for 
the 2;-aligned and planar x-aligned states of the current model in their Neel and stripe phases 
in the various LSUBn approximations are shown in Table [H Parallel computing is employed 
to solve the corresponding coupled sets of CCM bra- and ket-state equations (l6alb).— Our 
computing power is such that we can obtain LSUBn results for n = {2, 4, 6, 8, 10} for both 
the 2;-aligned model states and the x-aligned model states, as shown in Table [B However 
the very large numbers of fundamental configurations retained in the latter case at the 
LSUBIO level is only possible with supercomputing resources. For example, the solution of 
the equations involving the nearly 150,000 fundamental configurations for the stripe phase 
of the planar x-aligned state required the simultaneous use of 600 processors running for 
approximately 6 hours, for each value of the anisotropy parameter A in the Hamiltonian of 
Eq. (P. 
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TABLE I: Numbers of fundamental configurations (ft f-c.) retained in the CCM LSUBn approxi- 
mation for the z-aHgned states and the planar x-aligned states of the s = 1/2 Ji-^^-J2^^ model. 

z-aligned states planar x-aligned states 



Scheme (t f-C- tt f-c 



Neel stripe Neel stripe 



LSUB2 


1 


1 


1 


2 


LSUB4 


7 


9 


10 


18 


LSUB6 


75 


106 


131 


252 


LSUB8 


1287 


1922 


2793 


5532 



LSUBIO 29605 45825 74206 148127 



The final step in any CCM calculation is then to extrapolate the approximate LSUBn 
results to the exact, n — > oo, limit. Although no fundamental theory is known on how 
the LSUBn data for such physical quantities as the gs energy per spin, E/N, and the gs 
staggered magnetization, M, scale with n in the n oo, limit, we have a great deal of 
experience in doing so from previous calculations.— '^Si'^'^'^'^'^'^'^ Thus, we employ here 
the same well-tested LSUBn scaling laws as we have used, for example, for the Ji-J[-J2 
model,— 1^ namely 

E/N = ao + aiTi'^ + a^n'^ (10) 

for the gs energy per spin, and 

M = 6o + ^"°-^(&i + M"^) (11) 

for the gs staggered magnetization, both of which have been successfully used previously 
for systems showing an order- disorder quantum phase transition. An alternative leading 
power-law extrapolation scheme for the order parameter, 

M = Co + Cin-^2 , (12) 
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has also been successfully used previously to determine the phase transition points. For 
most systems with order-disorder transitions the two extrapolation schemes of Eqs. ffTTj) 
and fll2l) give remarkably similar results almost everywhere, as demonstrated explicitly, for 
example, for the case of quasi-one-dimensional quantum Heisenberg antiferromagnets with 
a weak interchain coupling.-^ However, in regions very near quantum triple points the form 
of Eq. (fTT!) is more robust than that of Eq. (fT2|) due to the addition of the next-to-leading 
correction term, as has been explained in detail elsewhere.— Hence, in this work we use the 
extrapolation schemes of Eqs. (ITOl) and (ITT!) . 

Obviously, better results are obtained from the LSUBn extrapolation schemes if the data 
with the lowest n values are not used in the fits. However, a robust and stable fit to any fitting 
formula with m unknown parameters is generally only obtained by using at least (m + 1) 
data points. In particular, a fit to only m data points should be avoided whenever possible. 
In our case both fitting schemes in Eqs. (fTOj) and (ITT]) have m = 3 unknown parameters 
to be determined. For all four model states we have LSUBn data with n = {2,4,6,8, 10}, 
and it is clear that the optimal fits should be obtained using the sets n = {4,6,8,10}. 
All the extrapolated results that we present below in Sec. IIVI are obtained in precisely this 
way. However, we have also extrapolated E/N and M using the sets n = {2, 4, 6, 8, 10} and 
n = {4,6,8}. In almost all cases they lead to very similar results, which adds credence to 
the stability of our numerical results and to the validity of our conclusions presented below. 

IV. RESULTS 

Figure [2] shows the extrapolated results for the gs energy per spin as a function of J2 (with 
Ji = 1) for various values of A, for the 2;-aligned and planar a;-aligned model states. For each 
model state, two sets of curves are shown, one (for smaller values of J2) using the Neel state, 
and the other (for larger values of J2) using the stripe state. As we have discussed in detail 
elsewhere,— "^"^ the coupled sets of LSUBn equations fl6al) have natural termination points 
(at least for values n > 2) for some critical value of a control parameter (here the anisotropy, 
A), beyond which no real solutions to the equations exist. The extrapolation of such LSUBn 
termination points for fixed values of A to the n, ^ 00 limit can sometimes be used as a 
method to calculate the physical phase boundary for the phase with ordering described by 
the COM model state being used. However, since other methods exist to define the phase 
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FIG. 2: Extrapolated CCM LSUBn results using the z-aligned and planar x-aligned states for the 
gs energy per spin, E/N, for the Neel and stripe phases of the s = 1/2 Ji-^^-J2'^^ model. The 
LSUBn results are extrapolated in the limit n — > oo using the sets n = {4, 6, 8, 10} for both the 
z-aligned states and the planar x-aligned states. The NN exchange coupling Ji = 1. The meaning 
of the E'max points shown is described in the text. 



transition points, which are usually more precise and more robust for extrapolation (as we 
discuss below), we have not attempted such an analysis here. 

Instead, in Fig. [21 the -Emax points shown, for each set of calculations based on one of the 
four CCM model states used, are either those natural termination points described above 
for the highest (LSUBIO) level of approximation we have implemented, or the points where 
the gs energy becomes a maximum should the latter occur first (i.e., as one approaches the 
termination point). The advantage of this usage of the -Emax points is that we do not then 
display gs energy data in any appreciable regimes where LSUBn calculations with very large 
values of n (higher than can feasibly be implemented) would not have solutions, by dint of 
having terminated already. 



Curves such as those shown in Fig. 2(a) illustrate very clearly that the corresponding pairs 
of gs energy curves for the 2;-aligned Neel and stripe phases cross one another for all values 
of A above some critical value, A > 2.1. The crossings occur with a clear discontinuity in 
slope, as is completely characteristic of a first-order phase transition, exactly as observed in 
the classical (i.e., s oo) case. Furthermore, the direct first-order phase transition between 
the z-aligned Neel and stripe phases that is thereby indicated for all values of A > 2.1, 
occurs (for all such values of A) very close to the classical phase boundary J2 = |, the point 
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FIG. 3: Extrapolated CCM LSUBn results using the z-aligned and planar x-aligned states for the 
gs staggered magnetization, M, for the Neel and stripe phases of the s = 1/2 Ji-^^-J2'^^ model. 
The LSUBn results are extrapolated in the limit n — > oo using the sets n = {4, 6, 8, 10} for both 
the z-aligned states and the planar x-aligned states. The NN exchange coupling Ji = 1. 



of maximum (classical) frustration. Conversely, curves such as those shown in Fig. 2(a) for 
values of A in the range 1 < A < 2.1 also illustrate clearly that the corresponding pairs of 
gs energy curves for the 2;-aligned Neel and stripe phases do not intersect one another. In 
this regime we thus have clear preliminary evidence for the opening up of an intermediate 



phase between the Neel and stripe phases. The corresponding curves in Fig. 2(b) for values 
of A < 1 tell a similar story, with an intermediate phase similarly indicated to exist between 
the xj/-planar- aligned Neel and stripe phases for values of A in the range —0.1 < A < 1. 

We show in Fig. [3] corresponding indicative sets of CCM results, based on the same four 
model states, for the gs order parameter (viz., the staggered magnetization), to those shown 
in Fig. [2] for the gs energy. The staggered magnetization data completely reinforce the phase 
structure of the model as deduced above from the gs energy data. Thus, let us now denote 
by Mc the quantum phase transition point deduced from curves such as those shown in 
Fig. El where is defined to be either (a) the point where corresponding pairs of CCM 
staggered magnetization curves (for the same value of A), based on the Neel and stripe 
model states, intersect one another if they do so at a physical value M > 0; or (b) if they 
do not so intersect at a value M > 0, the two points where the corresponding values of the 
staggered magnetization go to zero. 

Clearly, case (a) here corresponds to a direct phase transition between the Neel and stripe 
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phases, which will generally be first-order if the intersection point has a value M 7^ (and, 
exceptionally, second-order, if the crossing occurs exactly at M = 0). On the other hand, 
case (b) corresponds to the situation where the points where the LRO vanishes for both 
quasiclassical (i.e., Neel-ordered and stripe-ordered) phases are, at least naively, indicative 
of a second-order phase transition from each of these phases to some unknown intermedi- 
ate magnetically-disordered phase. We return to a discussion of the actual order of such 
transitions in Sec. |Vl In summary, we hence define the staggered magnetization criterion 
for a quantum critical point as the point where there is an indication of a phase transi- 
tion between the two states by their order parameters becoming equal, or where the order 
parameter vanishes, whichever occurs first. A detailed discussion of this order parameter 
criterion and its relation to the stricter energy crossing criterion may be found elsewhere.— 



From curves such as those shown in Fig. 3(a) we see that for A < 1.95 for the z-aligned 
states, there exists an intermediate region between the critical points at which M — for 
the Neel and stripe phases. Conversely, for A > 1.95 the two curves for the order parameters 
M of the quantum Neel and stripe phases for the same value of A meet at a finite value. 



M > 0, as is typical of a first-order transition. Similarly, Fig. 3(b) shows that for the 
planar x-aligned states, there exists an intermediate region between the critical points at 
which M — i> for the Neel and stripe phases for all values of A in the range —0.15 < A < 1. 
Again, the two curves for the order parameters M of the Neel and stripe phases for the same 
value of A intersect at a value M > for A < —0.15. In order to show more explicitly how 
the quantum phase transitions are driven by anisotropy. A, we display the same data for the 
extrapolated results for the order parameter, M, somewhat differently in Fig. HI where we 
plot M as a function of A for various values of J2 around the value J2 = 0.5, corresponding 
to the point of maximum (classical) frustration. 

By putting together data of the sort shown in Figs. [21 [3], and H] we are able to deduce 
the gs phase diagram of our 2D spin-1/2 J^^^-J^-^^ model on the square lattice, from our 
CCM calculations based on the four model states with quasiclassical antiferromagnetic LRO 
(viz., the Neel and stripe states for both the ^-aligned and planar xy-aligned cases). We 
show in Fig. [5] the zero-temperature gs phase diagram, as deduced from the order parameter 
criterion, and using our extrapolated LSUBn data sets with n = {4, 6, 8, 10}, shown as the 
critical value J| for the NNN exchange coupling J2 as a function of anisotropy A (with NN 
exchange coupling strength Ji = 1). Very similar results are obtained from using the energy 
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FIG. 4: Extrapolated COM LSUBre results using the z-aligned and planar x-aligned states for the 
staggered magnetization versus the anisotropy A for the s = 1/2 J^^^-J^^^ model, for the NN 
exchange coupling Ji = 1. The LSUBn results are extrapolated in the limit n — > oo using the sets 
n = {4, 6, 8, 10} for both the z-aligned model states and the planar x-aligned model states. 



criterion, where it can be applied (viz., along the transition lines between quasiclassical 
states with magnetic LRO). In order to test the accuracy of our results, particularly the 
positions of the phase boundaries shown in Fig. [5l we have also performed extrapolations 
using the LSUBn data sets with n = {2, 4, 6, 8, 10} and n = {6, 8, 10} for both the energy 
criterion and the order parameter criterion. In general terms we find that the results are 
remarkably robust, and the error bars quoted below are based on such an analysis. 

For the case of the z-aligned states, all of our results provide clear and consistent evidence 
for an upper quantum triple point (QTP) at (A'^ = 2.05 ± 0.15, J| = 0.530 ± 0.015) (for 
Ji = 1). For 1 < A < 2.0, there exists an intermediate paramagnetic (magnetically- 
disordered) quantum phase, separating the Neel and stripe phases. This intermediate phase 
disappears for A > 2.0, and both our energy and order parameter criteria give clear and 
unequivocal evidence for a direct first-order quantum phase transition between the two 
quasiclassical antiferromagnetic states in this regime, just as in the corresponding classical 
model (i.e., with s oo). The phase boundary approaches the classical line J| = 0.5 as 
A oo. 

Similarly, for the case of the xy-planar-aligned phases, a second (lower) QTP occurs at 
(A'^ = —0.10 ±0.15, J2 = 0.505 ±0.015) (for Ji = 1), with an intermediate disordered phase 
existing in the region —0.1 < A < 1. The z-aligned and xy-planar-aligned phases meet 
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FIG. 5: Extrapolated COM LSUBre results using the z-aligned and planar x-aligned states for the 
ground-state phase diagram of the s = 1/2 J^^^ -J^^^ model, for the NN exchange coupling 
Ji = 1. The LSUBn results for the staggered magnetization are extrapolated to the limit n — > oo 
using the sets n = {4, 6, 8, 10} for both the z-aligned model states and the planar x-aligned model 
states. Mc = magnetization critical point, defined in the text. 



precisely at A = 1, just as in the classical case. Exactly at the isotropic point A = 1, where 
the model becomes just the original J1-J2 model, the disordered phase exists for the largest 
range of values of J2, < J2 < can be clearly seen from Fig. [5l For the pure J1-J2 

model our calculations yield the values J2^/Ji = 0.44 ± 0.01 and J2^/Ji = 0.59 ± 0.01 that 
demarcate the phase boundaries for the disordered phase, in complete agreement with both 
our own earlier work and that of others that we have already discussed in Sec. [H 



V. DISCUSSION AND CONCLUSIONS 



We have shown in detail how, as expected, the quantum fluctuations present in the spin- 
1/2 J1-J2 model on the 2D square lattice, that has become an archetypal model for studying 
the interplay between quantum fluctuations and frustration, can be tuned by the introduc- 
tion of spin anisotropy. We have clearly confirmed our prior expectation that anisotropy 
reduces the quantum fluctuations. Thus, for both the cases A > 1 and < A < 1, the 
intermediate paramagnetic phase present in the pure J1-J2 niodel is observed to shrink to a 
smaller range of values of J2/J1 centered near to the point of maximal classical frustration, 
J2/J1 = |, that marks the classical phase boundary between the Neel-ordered and coUinear 
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stripe-ordered phases. 

We have seen that the intermediate disordered phase disappears precisely at two quantum 
triple points at = —0.10 ± 0.15 and A'^ = 2.05 ± 0.15, and that for values of A outside 
the range spanned by these values the intermediate phase is totally absent. In particular, 
for A > 2.0 we find unequivocal evidence for a first-order phase transition between the 
Neel and coUinear stripe phases. This direct first-order phase transition between states of 
different quasiclassical antiferromagnetic ordering is very similar to what has been observed 
in another similar extension of the spin-1/2 J1-J2 model on a square lattice, namely the 
so-called Ji-J2-J± model on a stacked square lattice where we now introduce a (weak) 
interlayer coupling through NN bonds of strength J±. The quantum fluctuations in the 
J1-J2 model are tuned here by the parameter J±. An analysis of this model^^ found that 
the intermediate region of disordered paramagnetic phase, a^^^ < J2/J1 < a^^, in the pure 
J1-J2 model now shrinks as the interlayer coupling strength is increased. The second- 
order phase transition for the Neel-ordered phase to the paramagnetic phase disappears for 
J±/Ji above some critical value (estimated to be in the range 0.2-0.3) marking a QTP in 
the J2-J± plane (with Ji = 1). Above the QTP there is again a direct first-order phase 
transition between the two phases of different quasiclassical antiferromagnetic LRO. 

On the other hand this scenario of a first-order phase transition between the two states 
of different quasiclassical LRO may be contrasted with the situation observed in yet another 
generalization of the pure spin-1/2 J1-J2 model on a square lattice, namely the so-called Ji- 
J[-J2 model that we have briefly discussed in Sec. [H In this case the quantum fluctuations 
are tuned by introducing a spatial anisotropy so that the NN bonds have different strengths 
in the intrachain (Ji) and interchain {J[) directions on the square lattice. A similar CCM 
analysis of the spin-1/2 version of this model by some of the present authors^^ again found a 
QTP in the J2-J1 plane (with Ji = 1), now below which the disordered paramagnetic phase 
disappears, and there is again a direct phase transition between the quasiclassical Neel and 
stripe-ordered phases with magnetic LRO. However, the surprising and novel situation found 
here was the existence of strong evidence for the phase transition in this case to be second- 
order, and hence inexplicable by standard Ginzburg-Landau theory, as discussed more fully 
in Sec. |T] above. 

Having discussed the transition line between the two phases of quasiclassical antiferro- 
magnetic LRO in the phase diagram in the J2-A plane of our spin-1/2 J^^^-J^^^ model 
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on the 2D square lattice, we turn our attention to the four phase boundary hues shown 
in Fig. [5] that dehmit the region of existence for the intermediate disordered paramagnetic 
phase. As has been explained in Ref. 59|| a judicious combination of the CCM energy data 
with the CCM order parameter data can shed light on the nature of the phase transitions 
between the quasiclassically long-range-ordered phases and the paramagnetic phase. The 
method for so doing relies essentially on the fact that although we perform our CCM cal- 
culations with model (or reference) states with quasi classical LRO, one knows^i^i^i^ that 
one can also reliably use such calculations in parameter regimes where all semblance of the 
quasiclassical LRO is destroyed. Thus, what is required for the CCM equations to converge 
to a solution is a sufficient overlap between the wave functions of the model (reference) state 
1$) and the true GS |\1'). The termination points of the CCM LSUBn equations discussed 
above are indicators of where this condition breaks down. Thus, provided that the CCM 
LSUBn equations converge and yield extrapolated solutions far enough beyond the points 
Mc where the order parameter vanishes, we can also determine whether the solution based 
on the Neel-ordered or the stripe-ordered model state has lower energy. 

We find in this way that there are indicators of a very narrow region where the gs energy 
obtained with the Neel model state might be slightly lower in energy than that obtained 
with the collinear striped model state, even in regions (close to) where the Neel order pa- 
rameter has already gone to zero, but where the stripe order parameter is still nonzero. 
As explained in more detail in Ref. 59|], the use of this evidence here points towards the 
zero-temperature phase transitions from Neel LRO to quantum paramagnetic disorder being 
second-order, while the transitions from quantum paramagnetic disorder to collinear stripe 
order are possibly (rather weakly) first-order rather than second-order. We stress, however, 
that the analysis here is very sensitive to the accuracy of our results, and the evidence for 
the nature of these quantum phase transitions involving the quantum paramagentic state in 
the regime —0.1 ^ A < 2.0 is less compelling than that for the transition between the two 
quasiclassically ordered states being first-order in the regime A > 2.0. 

The only other analysis of the current spin-1/2 J^^^-J^^^ model on the square lattice 
of which we are aware^ has been performed at the very low level of lowest-order spin-wave 
theory (LSWT). For the case studied here of equal spin-anisotropy parameters in the NN 
and NNN exchange bonds, these authors have only investigated the case A > 1, for which 
they find an (upper) QTP at a very small value of the anisotropy parameter, A^^^ ~ 1.048, 
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much smaller than the corresponding value A^^^ = 2.05 ± 0.15 obtained by us for the upper 
QTP. Such an extreme fragility or sensitivity of the paramagnetic phase to spin anisotropy 
is not easy to understand. In the face of our own much more accurate calculations it would 
seem simply to be an artefact of the LSWT approximation. On the other hand, the LSWT 
analysis does give the same qualitative trends as found by us for the phase transitions in the 
range A > 1, viz., a second-order transition between the Neel-ordered and disordered phases, 
and a first-order transition between the disordered and coUinear stripe-ordered phases for 
1 < A < A^^-j, and a direct first-order transition between the Neel-ordered and collinear 
stripe-ordered phases for A > A^^^. 

It is perhaps worth noting at this point in the context of spin-wave theory (SWT) that 
Igarashiii has shown that whereas its lowest-order (or linear) version (LSWT) works quite 
well when applied to the isotropic Heisenberg model with NN couplings only, it consistently 
oversestimates the quantum fluctuations in the pure (isotropic) J1-J2 model as the frustra- 
tion J2/J1 increases. Thus, he showed by going to higher orders in SWT in powers of 1/s, 
where LSWT is the leading order, that the expansion converges reasonably well for values of 
a = J2/J1 ^ 0.35, but for larger values of the frustration parameter a, including the point 
a = 0.5 of maximum classical frustration, the series loses stability. He showed for the s = ^ 
J1-J2 model that whereas LSWT predicts^ a value of a'^^ ~ 0.38 at which the transition 
from the Neel-ordered phase to the disordered phase occurs, the higher-order corrections to 
SWT for a < 0.4 make the Neel-ordered phase more stable than predicted by LSWT. This 
is precisely in agreement with our own predicted value of a'^^ = 0.44 ± 0.01 for the s = | 
J1-J2 model on the square lattice. He concludes that any predictions from SWT for the 
J1-J2 model on the square lattice are likely to be unreliable for values J2/J1 ^ 0.4. 



For reasons unclear to us, the authors of Ref. [48| never investigated the regime with 
A < 1, for which we find a lower QTP at A^^^ = -0.10 ±0.15, J^^^^ = 0.505 ±0.015. Clearly, 
our results are consistent with this lower QTP occurring exactly at the isotropic point 
(i.e., A = 0) of the model, and also exactly at the point of maximal classical frustration, 
J2 = |. A more detailed theoretical investigation of the corresponding J^^-J^^ model is 
clearly warranted by our results. 

Finally, we note that in our analysis here we have relied on two of the unique strengths of 
the CCM, namely its ability to deal with highly frustrated systems as readily as unfrustrated 
ones, and its use from the outset of infinite lattices. In turn, these features lead to its ability 
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to yield accurate phase boundaries even in the very dehcate regions near quantum triple 
points. Our own results for the gs energy and staggered magnetization from four sets of 
independent calculations based on different reference states provide us with a set of internal 

checks that lead us to believe that we now have a self-consistent and robust description of 
this rather challenging model system. 
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